DirectLiNGAMで構造方程式モデルを推定する
最初に
最近、統計的因果探索 (機械学習プロフェッショナルシリーズ) を読んでいて、データ生成過程における「構造方程式モデル」を把握する方法に興味を持ちました。
特に、本書籍の中で紹介されているLiNGAMについて調べていて、こちらのリポジトリで紹介されているDirectLiNGAMについて学んだので、備忘録も兼ねてブログにしようと思います。
また、本エントリーの内容は4.参照の内容を大いに参考にさせていただいているので、詳細が気になる方はそちらをご参照ください。
目次
1.はじめに
データ分析をしていると、よく「因果と相関は違うよね」という言葉を聞くと思います。
よく使われる説明として、国ごとの「チョコレート消費量」と「ノーベル賞受賞者数」が正の相関関係にあったとしてもそれだけでは因果関係まではわからない、というものがあります。
例えば、「チョコレート消費量」と「ノーベル賞受賞者数」が正の相関関係にある、といった時にそれらの関係は下記のように複数の因果関係が考えられますが、いずれの因果関係においても同じ相関係数が計算できる可能性があり、その意味するところは大きく異なります。
(他のパターンもあるかと思いますが、今回は省略します)
データ分析をするにあたっては、相関係数のみならず因果関係を知りたくなるケースが多いはずです。
そして、因果関係を探索するための1手法として「DirectLiNGAM」というアルゴリズムがあるので、それについて調べた内容についてまとめてみます。
2.DirectLiNGAMについて
DirectLiNGAMについてざっくり概説します。
ただ、その前に「構造方程式モデル」、「因果探索における3つのアプローチ」、「LiNGAM(Linear non-Gaussian acyclic model)」について概説します。
2-1.構造方程式モデル（Structural equation models: SEM）
「構造方程式モデル」は「変数の値の決定関係」を表すものです。
「データ生成過程を数式的に記述したモデル」と言ってもいいと思います。
これだけではよくわからないですね（笑）
「構造方程式モデル」を理解するために、簡単な例を踏まえてみます。
下記は「チョコレート消費量が増えると、ノーベル賞受賞者数が増える」といったパターンを表したものです。
先ほど見たものと比較すると、「e(x)」、「e(y)」という変数が増えています（外生変数）。
これらはそれぞれ、「チョコレート消費量(x)」、「ノーベル賞受賞者数(y)」を説明するために「必要だけど観測できなかった変数たちをまとめたもの（ex.地域/文化）」、と思ってください。
このケースは、数式的に言うと下記のように表現できます。
左辺は右辺から決定される（生成される）、という「データが生成される過程の構造がモデル化」できているわけです。
また、細かいですが構造方程式は下記4つ組から定義されます。
| 意味 | 上記の構造方程式で言うとどれ？ | |
|---|---|---|
| 内生変数 | 構造方程式モデルの左辺に登場する変数 | x、y |
| 外生変数（exogenous variable） | 構造方程式モデルの右辺にのみ登場する変数 | e(x)、e(y) |
| 内生変数と外生変数をつなぐ関数 | 内生変数を説明するために定義する関数 | f |
| 外生変数の確率分布 | そのまま | e(x)/e(y)を生成するところ |
なお、今回は明記していませんがe(x)/e(y)は「外生変数の確率分布」に基づいて計算されます。
なので、今回例示した構造方程式モデルでは、上記の4つ組み全てが登場しています。
2-2.因果探索における3つのアプローチ
さて、早速構造方程式モデルを求めていきたいところですが、その前に少しだけ考え方を整理します。
因果探索においては、「内生変数と外生変数をつなぐ関数」、「外生変数の確率分布」への仮定によって、構造方程式を求める際のアプローチ手法は「ノンパラメトリック」、「パラメトリック」、「セミパラメトリック」のいずれかに区分されます。
| 内生変数と外生変数をつなぐ関数 | 外生変数の確率分布 | |
|---|---|---|
| ノンパラメトリック | 仮定なし | 仮定なし |
| パラメトリック | 仮定あり | 仮定あり |
| セミパラメトリック | 仮定あり | 仮定なし |
いずれのアプローチにおいても、「因果的マルコフ条件（causal Markov condition）」や「忠実性（faithfulness）」を仮定します。
ノンパラメトリックなアプローチにおいては、一般的に「真の因果グラフ」を求めることはできませんが、「真の因果グラフと同じ条件付き独立性が成立する因果グラフの集合」を求めることができます。
また、パラメトリックなアプローチでも、「内生変数と外生変数をつなぐ関数が線形である」、「外生変数の確率分布がガウス分布である」という仮定をおいてしまうと、ノンパラメトリックなアプローチの際と同様に「真の因果グラフ」を一意に識別することができません。（ガウス分布の性質上）
2-3.LiNGAM(Linear non-Gaussian acyclic model)
LiNGAM(Linear non-Gaussian acyclic model)は、「内生変数と外生変数をつなぐ関数は線形」という仮定をおく一方、「外生変数の確率分布は非ガウス連続分布でそれぞれが独立（未観測共通原因なし）」という緩い仮定なのでセミパラメトリックアプローチとして分類されています。
具体的には下記のような式で表され、bは「パス係数行列」で、各変数から各変数への影響を保持した行列です。
（bは次元が「変数の数×変数の数」の正方行列で、対角成分が全て0である下三角行列です）
参照:構造方程式モデルによる因果推論:因果構造探索に関する最近の発展
xiは観測変数、eiは外生変数を意味しています。
この「パス係数行列」の「ゼロ/非ゼロパターン」を一意に識別できれば因果グラフを一意に推測できた、ということになります。
そして、LiNGAMモデルでは「外生変数の確率分布」の非ガウス性を利用することで、因果グラフを一意に推測可能です。
さて、LiNGAMモデルを採用する前に、「LiNGAMモデルの仮定を満たしているか」の検証を事前にすることもできます。(外生変数の非ガウス性、外生変数の独立性）
また、LiNGAMモデルを「未観測共通原因がある場合」に拡張したモデルに関する研究も進んでいます。
本エントリーでは対象外としますが、「未観測共通原因がある場合」の方が実際には多いはずなので、そちらも今後ブログ化していこうと思います。
2-4.DirectLiNGAMについて
「DirectLiNGAM」は「変数間の順序関係」がわかるまで「変数間の単回帰」、「独立性評価」を繰り返すことでLiNGAMモデルを推測しよう、というものです。
逐次、「もっとも外生変数らしい」変数の影響を他の変数から削除していくことでパス係数行列を求める、というものですね。
こちらのp.95~p.104がとても参考になります。
また、もし「事前知識(ex.どの変数がどの変数に影響があるか)があって、それを適用したい」という場合は事前知識をモデルに反映することで「より正確な結果を求めつつ、計算量を減らすこと」も可能です。
詳細はDirectLiNGAM: A direct method for learning a linear non-Gaussian structural equation modelの「3.4 Use of prior knowledge」をご確認ください。
3.やってみる
前置きが長くなってしまいましたが、実際にこちらを参考にDirectLiNGAMを試してみようと思います。
import numpy as np import pandas as pd import graphviz import lingam from lingam.utils import make_dot print([np.__version__, pd.__version__, graphviz.__version__, lingam.__version__]) np.set_printoptions(precision=3, suppress=True) np.random.seed(100)
['1.19.5', '1.1.5', '0.16', '1.5.1']
今回使うライブラリのバージョンはこんな感じです。
x3 = np.random.uniform(size=1000) # 一様分布 x0 = 3.0*x3 + np.random.uniform(size=1000) # 一様分布 x2 = 6.0*x3 + np.random.exponential(2,1000) # 指数分布 x1 = 3.0*x0 + 2.0*x2 + np.random.gamma(2., 2., 1000) # ガンマ分布 x5 = 4.0*x0 + np.random.uniform(size=1000) # 一様分布 x4 = 8.0*x0 - 1.0*x2 + np.random.beta(2, 2, 1000) # ベータ分布 X = pd.DataFrame(np.array([x0, x1, x2, x3, x4, x5]).T ,columns=['x0', 'x1', 'x2', 'x3', 'x4', 'x5']) X.head()
ほぼほぼサンプルをそのまま利用しますが、そのままサンプルをやっても微妙なので、少し外生変数の確率分布を弄ってみます。
さて、こんな感じのデータに対してDirectLiNGAMで構造方程式モデルを推測してみます。
model = lingam.DirectLiNGAM() model.fit(X) make_dot(model.adjacency_matrix_)
いい感じですね。
